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1. Introduction 

The interaction of optical fields with Bose-Einstein condensates (BECs) has attracted 
considerable interest over the last years. A rather intriguing phenomenon in this context 
is the generation of superradiant light pulses accompanied by condensate momentum 
side-modes, when far off-resonant laser light is scattered from an elongated BEC 
[Il[2l[3llll[5l[6l[71[8]. The process is initially driven by quantum fluctuations in 
the vacuum electromagnetic field (i.e., spontaneous Rayleigh scattering) [6]. After the 
first photons have been scattered, bosonic stimulation sets in and additional successive 
scattering events give rise to a rapid growth of the emitted radiation and of the side-mode 
populations. Moreover, due to the geometry of the condensate, the gain is largest when 
the scattered photons leave the condensate along its long axis, in the so called endfire 
modes. The fact that the process is triggered by quantum fluctuations is expected to be 
reflected in various macroscopic parameters of the superradiant (SR) light pulses (e.g., 
peak intensity, shape, etc), which should also fluctuate when the same experiment is 
repeated. 

Throughout this work we focus on the passage-time statistics, which refer to the 
fluctuations in the time at which the pulse intensity at a given point in space reaches a 
certain reference intensity. The effect of quantum fluctuations is expected to be rather 
prominent for short times, where the observed intensity is relatively small. Given that 
the populations of the endfire modes and the side-modes grow fast in time, the evolution 
of the system becomes practically deterministic very soon |8] . For our purposes, it 
therefore suffices to focus on the early stage of the phenomenon, where the BEC is 
practically undepleted, only first-order atomic side-modes are significantly occupied, 
and the dynamics of the system are governed by linearized equations of motion. In view 
of the above arguments, the passage-time statistics at the early stage are expected to 
represent adequately also the delay-time statistics, which refer to the first maximum of 
the pulse that appears at later times. 

Analogous theoretical as well as experimental studies have been performed in the 
context of conventional (Dicke) superradiance [9l[T0l[Tl], where excited atoms interacting 
with the vacuum of the electromagnetic field emit collectively; analytic expressions 
for the probability distribution of the passage and delay times have been obtained 
p!T| [12] . In the system under consideration, however, condensed atoms may also scatter 
a photon from the endfire modes back into the laser mode, giving thus rise to recoiling 
atoms moving against the direction of the applied laser field (backward atomic side- 
modes). This essential feature of superradiance from Bose-Einstein condensates is a 
key difference between this process and conventional superradiance. In this work we 
investigate how the presence of backward atomic side-modes affects the passage-time 
statistics, and examine whether the corresponding probability distributions usually 
discussed in conventional superradiance can also describe accurately the fluctuations 
in systems involving ultracold atoms. To the best of our knowledge, none of these 
questions has been addressed either theoretically or experimentally so far . In the 



Passage-time statistics of superradiant light pulses from Bose-Einstein condensates 3 



following section we present the theoretical framework, while in Sec. [3] we analyze the 
passage-time statistics and their scaling behaviour. We conclude in Sec. HI 

2. Theoretical Framework 

The system under consideration pertains to a BEC elongated along the ^-axis, consisting 
of atoms. A linearly polarized laser with frequency ui = kic, far detuned from the 
closest atomic transition by a value of 6, is illuminating the cloud along the x-axis. 
Due to the coherent nature of the condensate, successive Rayleigh scattering events 
are strongly correlated and lead to collective superradiant behaviour. As a result of 
the cigar-shape of the condensate, the gain is largest when the scattered photons leave 
the condensate along its long axis, traveling in the so called endfire modes with wave 
vectors k i/c^e^. A condensate atom can scatter a laser photon into the endfire 
modes, experiencing a recoil ~ h{kiex — k). On the other hand, it can also scatter a 
photon from the endfire modes into the laser mode, in which case its momentum changes 
by ~ h{—kiex + k). These processes lead to the formation of two pairs of atomic 
side-modes, consisting of counterpropagating atoms with a narrow momentum spread 
(compared to ki). Of course, atoms within these side-modes can also scatter photons, 
thereby acquiring higher momenta, but since we are interested in the early stage of the 
process, we consider only first order side-modes to be populated. 

2.1. The Model 

Neglecting any coupling between counterpropagating photonic endfire modes, the 
dynamics in the two halves of the elongated system become independent . We can thus 
focus on the endfire modes with k -\-kez and consequently on the two sets of atomic 
side-modes (with central momenta ^q = ±h{kiex — ke^)) that are coupled to them. Due 
to the strong confinement along the x- and y-axis, we can assume the transverse profiles 
of the matter field ip±{x, y) to be well described by a classical function, independent of 
the 2;-coordinate [7]. Assuming a Fresnel number close to unity for the electromagnetic 
fields, we can apply the same approximation for the transverse part of the radiation 
field u±{x,y), effectively reducing the problem to one dimension. 
We expand the field operators as 

1 

^(x,t) = M^,y) J2 Uz,t)e'^^''^-'^'>-'^^' (la) 
i=-i 

EW(x,t) = ^eye'('^'^-^''^+u^{x,y)Ei^\z,t)ey, (16) 

where u±i = h{kf -\- k^)/2M with M the atomic mass and uq = 0. The matter- wave 
operator is split up in three parts, describing the two side- modes (j = ±1) and the BEC 
at rest (j = 0), 

i:j{z, t) = e'-^* -7fC_,k+p{t), (2a) 
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where Aq is the interval {—k/2,k/2) in /c-space, L is the length of the BEC and Cp 
annihilates an atom with momentum hp. Since the BEC at rest remains practically 
undepleted it can be treated as a time independent classical function and hence we can 
set ipoi^yt) = i^oi^)- Similarly, we expand the endfire mode operator as 



E 



2eo ^fl 




{2h) 



where Op is the photon annihilation operator. The frequencies Uk+p are approximated 
by a; = /c/c, the frequency of the scattered photons, and can therefore be taken out of 
the sum. This approximation is justified by the fact that dominant contributions to the 
sum come from momenta of order which is several orders of magnitude smaller 
than k. 

The one-dimensional field operators satisfy usual boson commutation relations 



e+{zi,t),e\{z2,t) 
e+{z,ti),e\{z,t2) 



5ij5i^{zi - Z2) 
SaIzi - Z2), 



1. / 
c 



t. 



(3 a) 

m 

(3c) 



where 6ij denotes the Kronecker delta and S^iz) is a delta-like distribution with width 
of order 1/k. Since the fields vary slowly on length scales of the order of 1/k, we 
will approximate the Sa{z) by actual delta functions. Inserting expansions (12. ip in 
the Maxwell-Schrodinger equations for the coupled matter-wave and electric fields and 
applying the slowly-varying-envelope approximation, we find the equations of motion 
for the operators 

5^li(e,r) 



dr 

dr 

<9e+(^,r) 5e+(^,r) 

+ x- 



dr 



iKe+(^,r)?/^o(0> 
-i«:e+(e,r)^o(0-2iV'-i(e,r) 



(4a) 
(4c) 



Here we have defined ip'_i{C,,T) 
dimensionless units 



e ijj-i{^,T) and rescaled length and time to 



T = 2uJrt, 



C = kiz, 

where Ur = hkf /2M. Accordingly, the fields are rescaled as 

e+(^,r) = -^e+{z,t) V^j(^,r) = -^iijj{z,t), 



(5) 
(6) 
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and the speed of light becomes x = The effective ID-couphng is given by 

K = g^/kiL/ {2uJr) with 



For a detailed derivation of fl2.ll) see [S]. 

Backwards recoiling atoms are a particular feature of superradiant Rayleigh 
scattering off condensates. The physical process underlying the backwards modes 
violates energy conservation by an amount AE ~ 4^^- Thus, according to Heisenberg 
uncertainty principle, it can take place only for times shorter than a critical time h/ AE 
which, in our units, is given by P, 

r, = 0.5. (8) 

For such short pulses, and for sufficiently high power one typically observes an X-shaped 
pattern for the distribution of the atomic side-modes with the initial BEC in the center 
and the recoiling atoms moving both in and against the direction of the applied laser 
pulse [Kapitza-Dirac or strong-pulse regime). On the other hand, for weaker pulses 
with duration longer than Tc, the distribution of the side-modes exhibits a fan pattern, 
involving mainly forward recoiling atoms [Bragg or weak pulse regime). If we neglect the 
atomic backwards side-mode altogether, (12. ip become formally equivalent to descriptions 
of "conventional" superradiance in excited atomic gases [9]. 

2.2. Solutions to the equations of motion 

We can use the Laplace transform to find exact solutions to the system (12. ip in terms of 
the operators evaluated at the boundary of their domain - i.e. at ^ = and r > and 
at ,^ > and r = 0. More details on this procedure are given in [8]. In the following we 
are interested in the radiation field only, whose time evolution is given by 



e+(e,r)= / dr'e+(0,r')Fo,o(75,o,r-r') 
Jo 

-- l^de k(n^li(e',o)Fi,o(7a',r) 

X Jo 



X Jo 

+Vo(e')^-i(e',0)Fo,i(75,e,T)] , (9) 
where we have neglected retardation effects and introduced 

ii,e = ^m-piai (10) 

X 

with p(^) = d^'l^/^ol^') The functions F^^y{u,v) are defined as 

where denotes the inverse Laplace transform. Explicit expressions for the functions 
F^^u{u,v) appearing in ([9]) are given in the appendix. They are combinations of 
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Bessel functions and integrals thereof. For the numerical calculations, we assumed 
the BEC to consist of = 10^ Thomas-Fermi distributed *^Rb atoms, such that 
^0 = ^/e{z)N6{Lz - z^)/L^ with L = 130/im. For the incoming laser we chose 
a rectangular profile and a wavenumber ki = 8.05 x 10^m~^, which results in a 
dimensionless length of the BEC of A ~ 1000. Coupling strengths are conveniently 
expressed in terms of the superradiant gain T = hi^N/ x, whose value separates the two 
regimes identified by experimental observations of superradiance from condensates [3]. 
Typically, the weak coupling regime is characterized by ~ lO^s^^ and F ~ 1, while 
g ~ lO^s"^ and F ^ 1 in the strong coupling regime. In our calculations, we chose 
F = 1 and F = 100 for the two regimes. 

2.3. Averaging over semiclassical trajectories 

The presented model is closely related to the mean- field (MF) model presented in [7]. 
In fact, we arrive at the MF equations with the same approximations by replacing all 
operators in (12.11) by their expectation value. Due to the generality of the Laplace 
transform, the same replacement in ([9]) gives the solutions to the MF model. A detailed 
comparison between the two models is given in [8]. One of the problems of the MF 
model is the fact that it cannot describe the spontaneous onset of the superradiance 
process. In order to start the evolution, one has to "seed" one of the modes with some 
initial values. In the context of conventional superradiance, Haake et al. introduced the 
idea of averaging over many semiclassical "trajectories" to obtain quantum results |12j . 
The MF equations are initially seeded with random variables according to a particular 
distribution. To obtain a particular quantum expectation value, one has to average over 
various solutions for the corresponding semiclassical quantity. We have investigated the 
extension of this idea to superradiant Rayleigh scattering off BECs where, in contrast 
to conventional superradiance, backwards recoiling atoms are also present. 

To see how this works, let us assume we want to calculate the normal-ordered nth 
order correlation function of the emitted light. 

[e+(e',rO]"). (12) 

From ([9]) and using the commutation relations as well as the fact that our initial state 
is the vacuum for all modes, we find the only non-vanishing expectation value involved 
to be 

(^+i(e(^\o) . . .^+i(e(^),o)^ii(e("+^\o) . (is) 

where the ^'•■'^ are integrated from to for j = 1,...,?7, and from to ^' for 
j = n + 1, . . . , 2n. Using the commutation relations, correlation ( fT3l) reads 

n 

TT j = l 

where 5 denotes the Dirac delta function and the sum runs over all permutations vr of 
order n. 



( eUe,r) 
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Figure 1. (Colour online) Example of the photon density within the BEC volume 
obtained through averaging over randomly seeded MF solutions. The quantum solution 
is shown in solid black, while averages over 20 and 2000 trajectories are shown in red 
dashed and blue dotted lines. Parameters: F = 1, and r = 2. 



Let us now seed the semiclassical model with [8] 

V^H_i(e,0) = Q, (15a) 
V'-i(^,0) = 0, (15&) 
e+(e,0) = 0, (15c) 

where is a random, normally distributed complex variable, with zero mean 
and variance 1 / a/A^, with the spatial step of a numerical implementatioiJ||. 
The average of the product ^+i(^(^\ 0) . . . V^+i(^("\ 0)ipl^{^^''+^\ 0) . . . ^;i(^(2"\ 0) over 
many trajectories (seedings) will effectively converge towards (IMI) . Due to the formal 
equivalence of the semiclassical solutions to the quantum ones, the average of the product 
e+(^i, T)"e+(^25 t)"" will consequently converge towards the quantum expectation value 
( 1T2|) . To find correlation functions of other operators, other seedings have to be used, 
which can be found in an analogous way and are summarized in [8]. 

As far as densities and populations are concerned, the convergence of the averaging 
procedure towards the quantum solution is fairly fast. In figure [H we plot the photon 
density 

n(^,r) = (eU^,^)e+(^,r)) (16) 

averaged over different numbers of trajectories compared to the quantum result. 
Typically, one has to average over a couple of thousand trajectories, to obtain well 
converged, smooth density profiles. Clearly, the averaging procedure is also applicable 
to correlations of higher order. One can calculate any normal ordered correlation 
function of the system by averaging over a sufficiently large ensemble of MF trajectories. 



I In our case 400 spatial gridpoints 
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provided that the operators involved in the correlation function have the same seeding 
requirements. The convergence, however, becomes rapidly slower with every order 
added, such that higher order correlations will require a larger number of trajectories. 
Higher order correlations are related to higher order statistical moments of a particular 
observable, such as the variance or skew of a particular distribution. 

3. Passage- Time Statistics 

We proceed to apply the technique of averaging over semiclassical trajectories to the 
study of passage-time statistics. First, we will find a simple probability density function 
which matches the distributions. Following this, we will investigate the behaviour of 
the distributions with respect to scaling of the two most important parameters of such 
an experiment and find invariant quantities for both cases. 

3.1. Distribution of Passage- Times 

In view of Sec. 12. 3[ a single semiclassical trajectory represents a single experiment on 
superradiance off a given BEC P, [TT], [12]. The process is initially driven by quantum 
fluctuations; within the MF model this is simulated by seeding ip^i with Gaussian white 
noise. After the first spontaneous scattering events, bosonic stimulation leads to a rapid 
growth of the emitted radiation and of the side-mode populations. The signature of 
the quantum fluctuations which started the process, however, is imprinted onto the 
superradiant light pulses. 

One way to investigate quantum fluctuations in this context, is by studying the first 
passage-time statistics, which pertain to the distribution of times at which the intensity 
of the SR pulse reaches a certain reference intensity. Equivalently, one may look at the 
times for which the number of photons that have left the condensate volume, given by 



reaches a certain threshold Performing a large number of semiclassical simulations, 
each one seeded with a new random tp+i, and following the evolution of ^(r) for 
every trajectory, one can record the corresponding passage-times T. The quantum 
fluctuations responsible for the start-up of the system are translated into fluctuations 
of the passage times, whose statistics can provide information about the onset of the 
phenomenon. Moreover, due to the rapid growth of the mode populations and the 
deterministic evolution of the system, the passage-time statistics are expected to also 
reflect the delay-time statistics, which refer to the times at which the first maximum of 
the SR pulse occurs. 

Based on these ideas, Haake et. al. investigated the statistics of passage- (and 
delay) times of SR light pulses obtained from a conventional superradiant system [T2] . 
A semi-analytic expression was given for the distribution of the passage-times, which 




(17) 
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was found to be well approximated by the Gumbel probability density function (PDF) 



Their predictions about the standard deviation of the passage-times were found to be 
in good agreement with experimental data |13] . 

In contrast to conventional superradiant systems, in superradiance off BECs we 
also have the presence of backwards atomic side-modes. This is an essential feature of 
our system and can influence statistical properties of the scattered light, such as the the 
first- and second-order correlation functions [8]. Here, we will discuss its influence on the 
passage-time statistics. To this end, we performed a number of numerical simulations 
using the semiclassical equations of motion. For each trajectory, the equations were 
seeded in accordance with (12. 3p . and we followed the time-evolution of the number of 
outgoing photons ^ij). As soon as the threshold was reached, we recorded the 
corresponding time t = T, and started a new simulation with new random seeding. The 
histograms depicted in figure [2] describe the distribution of the recorded times after a 
large number of simulations. 

In the weak pulse regime, the distributions are narrower (i.e., the variance of the 
distributions is smaller) than in the strong pulse regime. An intuitive understanding for 
this broadening can be obtained by recalling that in the strong pulse regime, the growth 
of the photonic endfire modes is slowed down by the presence of backwards recoiling 
atoms, as each backscattered atom removes a photon from the endfire modes. Therefore, 
with respect to the scaled time Fr, the spread of solutions for ^{j) corresponding 
to different seedings, is larger in the strong than in the weak pulse regime. This is 
illustrated in figure |3l where the evolution of ^ij) is plotted for two different seedings 
in the strong- and the weak-pulse regimes. 

All the distributions exhibit positive skew, which limits the PDFs appropriate to 
model the data. Several two-parameter, positively skewed PDFs were fitted to the data 
according to a least-square fit and assessed by calculating the residual errors as well as 
comparing mean, variance and skewness of the distributions to the ones of the data. 
While the Gumbel PDF of conventional superradiance gives a very satisfying fit, it is 
for all parameter values studied inferioi[§| to a fit with the inverse Gaussian PDF 



which, as a further advantage, does not support negative T. If a random variable X is 
distributed according to f[T^ . the mean and variance of X are given by 



In figure [2] we show two fits to the numerical data, corresponding to the Gumbel and 
the inverse Gaussian PDFs. 

§ The assessment of fitted PDFs was performed by calculating the sum of residual errors as well as by 
checking the agreement of the PDFs mean, variance, skew and kurtosis with the ones of the raw data. 



m 



PG(T;/i,A) = A"V^e-^P(-^). 





(19) 



X = /i, var(X) = ^?/\. 



(20) 
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Figure 2. (Colour online) Passage-time statistics of light pulses in superradiance from 
BECs. Histograms (a) and (c) show the distribution of the times at which the number 
of photons that have left the BEC reaches the threshold = 10, for (a) F = 1, and 
(c) r = 100. The corresponding histograms with neglected backwards recoiling atoms 
(equivalent to conventional superradiance), are plotted with inverted ordinate (b,d). 
Dark (red) histograms show the difference between the two cases, in the two regimes. 
The curves give least square fits with inverse Gaussian (solid) and Gumbel (dashed) 



We now turn to the comparison of the passage-time distributions for superradiance 
off BECs, to the ones obtained for conventional superradiance. As discussed in section 
12. H the latter system can be simulated within our model by neglecting the backwards 
atomic side- modes in (12. ip . The corresponding distributions of the recorded passage- 
times are depicted in figures [2] (b) and (d). From the difference of the two histograms, 
also plotted in figure [21 we see that in the weak pulse regime (T = 1), the two 
distributions are virtually the same [compare figures |2] (a) and (b)]. For weak couplings, 
the number of backwards recoiled atoms is not big enough to change dramatically 
the passage time statistics. In the strong pulse regime, however, the distribution for 
conventional superradiance is narrower than the one for superradiance off BECs, while 
it is also shifted towards smaller passage times [compare figures [2] (c) and (d)]. Finally, 
in the weak pulse regime, fits with either of the aforementioned two PDFs agree very 
well with the distributions of the numerical data [see figure |2] (a) and (b)]. On the 





PDFs. 
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Figure 3. (Colour online) Evolution of the number of photons in the SR light pulse 
for two different seeds, as a function of scaled time in the strong and weak coupling 
regime. The vertical lines indicate the wider spread due to the slower growth in the 
strong coupling regime. 

contrary, in the strong-pulse regime, and if backwards recoiled atoms are present, we 
find noticeable discrepancies between the two fits, particularly for small T. They become 
more prominent for small thresholds 

3.2. Scaling with coupling strength 

We have seen that our equations of motion can describe conventional superradiance if we 
set all the couplings to backward atomic side-modes to zero. The resulting reduced set 
of equations is also expected to be a good approximation in the weak-coupling regime 
of superradiance off condensates, provided that the time scales of interest are such 
that the backward atomic side-modes are scarcely populated. Assuming furthermore 
a homogeneous condensate, analytic solutions for the operators can be found, which 
moreover are a function of the scaled time Fr rather than just the time alone [8]. If 
one multiplies inverse- Gaussian distributed random numbers with a constant factor, 
the resulting numbers are again inverse- Gaussian distributed, with the mean /i and 
the shape parameter A rescaled by the same number. Hence the scaling behaviour is 
also reflected in the distributions, which can be seen by comparing the distributions of 
figure [2] (b) and (d), where a scaling factor or 100 is clear. We can therefore extend 
our previous observation on the scaling behaviour of the passage-time distributions for 
conventional superradiance, to the context of superradiance off condensates and define a 
more general "scaled" version of the inverse Gaussian PDF, which describes accurately 
passage-time distributions for various couplings within the weak-coupling regime, 

Vwp{r; /i. A) = Vig{t, /iwp/r, Awp/F). (21) 
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For any o^th, the two open parameters /iwp and Awp are constants characterizing the 
weak-pulse regime. For A/'th = 10, their value was estimated to /iwp ~ 3.2 and Awp ~ 16. 
In the case of conventional superradiance, it is known that the mean of the passage 
times scales as a constant divided by the atomic decay rate, while their variance scales 
as another constant divided by the decay rate squared [H]. In our model, the quantity 
corresponding to the decay rate is the superradiant gain F, which is directly proportional 
to the superradiant scattering rate. From fl20|) and fl2T]) we therefore find that the inverse 
Gaussian distribution recovers the same scaling laws as found in |llj . 

If we increase the coupling strength and leave the weak-coupling regime, we lose 
this scaling behaviour, as can be seen by comparing the histograms (a) and (c) in 
figure [2j Surprisingly, however, we recover fl2T|) deep within the strong pulse regime, 
albeit with new constant parameters, which for =yfth = 10 are given by Hsp ~ 5 and 
Asp ~ 20. To gain further insight into the change of the scaling behaviour of the 
passage-time distribution with the coupling strength, we calculated numerically the 
distributions ?^(FT, ^A^h), of the rescaled passage-times FT, for various values of F with 
Tiow < r < Thigh, and various thresholds o/ffh- The upper and lower values of F were 
chosen as Flow = 0.1 and Fhigh = 100, so that both regimes are covered. Subsequently, 
for a given we calculated the normalized overlap 



and the corresponding density plot is depicted in figure HI While the parameters /xwp/sp 
and Awp /sp are independent of the coupling, they are nontrivial functions of the photon 
threshold and their numerical dependence is also depicted in figures S] (a,c). Quite 
obviously, the scaling law within the two regimes is very well obeyed by the system. We 
can therefore conclude that in the early stages of superradiant scattering off condensates, 
the distribution of passage times in either regime is universally well modelled by the 
inverse Gaussian PDF and two constant parameters. The transition between the two 
scaling behaviours begins in the region, where the parameters and F match up 
in a way such that T ^ Tc- The distributions in this transient regime do not differ 
qualitatively from the ones presented in figure Ej but one cannot define a universal 
scaling law for them. As can be seen in figure H] (b), however, this transient regime (dark 
region) in the parameter space is rather small compared to the ones where the system 
is unambiguously found in either of the two regimes (white regions). As a final remark, 
it should be noted that for any threshold ^h? the ratio of mean and shape parameter 
of the resulting distribution is independent of F, which by virtue of ( l20i) implies that 
the ratio of the variance and the squared mean of passage times is independent of F. 

3. 3. Scaling with threshold photon number 

Having investigated the behaviour of the passage time distributions for varying F, it is 
natural to ask how they behave with a constant coupling but varying the other significant 
parameter of the system, the photon threshold ^4^h- For any choice of F, we find again a 




(22) 
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Figure 4. (Colour online) (a) Shape parameters of the inverse gaussian PDFs as a 
function of photon threshold yKh- (b) Density plot of the overlap p2|) . for various 
couplings r and thresholds White corresponds to a value of unity. Note the 

logarithmic scale for F. The contour lines indicate the points for which the expectation 
time for T is equal to 0.5,1 and 2 respectively. Each histogram contains 10'' trajectories, 
(c) Means /i of the scaled passage times as a function of photon threshold 



(23) 



simple combination of shape parameter and mean, which remains constant independent 
of ./({h- It is the ratio 

_ A _ f 
/x^ var(T) ' 

where the second equahty is a simple consequence of (l20l) . Figure [5] (a) depicts a 
contour plot of the value of s for different couplings and photon thresholds, and reveals 
the independence of this value of the coupling T. Since these results are based on the 
analysis of a large ensemble of random trajectories, a certain variation is to be expected, 
and the values of s as a function of F are obtained by averaging s over y({h- Moreover, 
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Figure 5. (Colour online) (a) Contour plot of the value of s as defined in ((23)) . (b) 
The parameter s, averaged over -yith, as a function of the coupling F (crosses) and the 
fits (24) in the weak and strong pulse regimes (lines). For every threshold and coupling 
f 0^ trajectories were evolved. 



for each regime, the quantity s has a very simple hnear dependence on the couphng F, 
which is illustrated in figure ([5]) (b). The functions plotted are the fits 

S'^p(r) = 1.56r (weak pulse regime) (24a) 

sgp(r) = 0.7ir + 2.86 (strong pulse regime). (246) 

Note that the equation for the strong pulse regimes is valid for large values of F, such 
that the offset becomes negligable and sgp is directly proportional to F. Hence this 
relation is consistent with equation[2Tl Moreover, also in this case the transition between 
the two fits (24), takes place in a regime where the parameters result in times around 
critical time Tc = 0.5. 

Random variables distributed according to the inverse Gaussian PDF and with 
this very ratio as a constant, also appear in different contexts. Consider the stochastic 
process 

Xo = (25 a) 

Xt=i^t + aBt, {25b) 

where Bf is the standard Brownian motion and > is a positive, constant drift. The 
times at which Xt crosses a certain limit a > are distributed according to the inverse 
Gaussian distribution [13] 

r„ = inf{0 < t\Xt = a} 



^IG T,;/i=-,A = — . (26) 



The ratio of shape parameter and squared mean is independent of a, which is analogous 
to s being independent of the photon threshold and we can thus make a connection 
between Brownian motion and superradiance off BECs. In particular, we can construct a 
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stochastic system analogous to (13. 3p . but in terms of T and s. Note that the particular 
choice of corresponding parameters is not unique. From f l26|l we find that a/u must 
correspond to the quantum expectation value T, while 

s = ^. (27) 

We can therefore set the drift u to unity and find T = a and a = 1/ ^/s. Hence, the 
times the stochastic system 



will cross Xt = T will reproduce the distribution of passage-times for a superradiant 
BEC emitting photons. This result by no means implies that superradiance 
is in fact a physical realization of Brownian motion. In contrast to Brownian 
motion, superradiance off condensates is a deterministic evolution from random initial 
conditions. After its initial stage, the system can be described by deterministic 
equations. This is never the case for Brownian motion. 

4. Conclusion 

We have used random seedings to investigate the passage-time statistics of superradiant 
light pulses emitted from an elongated atomic condensate. In general, the distributions 
of the passage times are very well approximated by the inverse Gaussian probability 
density function. The distributions in the Kapitza-Dirac regime, however, appear wider 
than in the Bragg regime, while they are also shifted towards larger passage times. 
This is due to the presence of backwards recoiling atoms in the system, which in the 
Kapitza-Dirac regime lead to a slow growth of end-fire modes. For both regimes the 
mean of the distributions is inversely proportional to the superradiant gain F, whereas 
their variances scale with F~^. The corresponding proportionality constants are different 
in the two regimes and for each regime, the ratio of shape parameter A and mean /i is 
independent of F. For a small region in parameter space, the system is not in either of 
the two regimes but finds itself in a transient one, where the distribution of the passage 
times changes smoothly between the two reference (outer most) distributions. This 
transient region is characterized by timescales around the critical time, at which energy 
conservation suppresses backwards recoiling atoms. For constant coupling and varying 
threshold, the ratio A//i^ is constant. The latter behaviour leads to a phenomenological 
correspondence between superradiance from BECs and Brownian motion with drift. 
We focused at the very early stage of the phenomenon, but due to the deterministic 
evolution of the system the observed behaviour can be expected to hold for later times 
and consequently the distributions are expected to also describe the delay-times of a 
superradiant pulse from a condensate. Finally we want to stress that the presented 
results are accessible to experimental verification. Similar experimental studies have 
been performed for conventional superradiance [TH [13] . 



^0 = 0, 

Xr = T+ -^B. 





(29) 
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Appendix A. Inverse Laplace Transforms 



Using elementary techniques for Laplace inversion, one finds 

F,^o{y,z) = io{2^)e{z) 

pz —Viz' 

Fo,i(|/,^) = e-2^Vo(2V^)e(^) 



-i2z' 



-A 



z — Z' 



F,4y,z) = e{z) / d^'e-2^7o 



,z — z 



2^y{z-z' 
Jo (2^/yz' 



Jo (2v^) (1.16) 



:i.lc) 



I z — z' 



-Q{z) dz'e-''^ J^^h[2^y{z-z')]M2^yz'] (l.ld) 



Fo,2{y,z) = e-''\l-Ji{2^)e{z) 



Q(z) / dz'e 



/~2iz' 



Z' 



-h[2^y{z-z')]Ji[2^] (Lie) 



z — z' 

with Ji and Jj the ith Bessel function of the first kind and the ith modified Bessel 
function respectively [8]. 
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